# coding:utf-8
from numpy.core.numeric import zeros, array, dot
from numpy.ma.core import sqrt, size, sin, cos, log10
from numpy.lib.function_base import append
from numpy.core.fromnumeric import transpose
from numpy.linalg.linalg import solve
import math, sys, pylab

class Cylshlfm(object):
    
    ### 境界条件 ###
    free = 0
    ssupport = 1
    gsupport = 2
    pointload = 3
    lineload = 4
    surfaceload = 5
    axisymln = 6
    axisymsf = 7
    psupport = 8
    ##############
    young = 1.e10 #ヤング率
    row = 1000.0 #密度
    poisson = 0.1 #ポアソン比
    thick = 0.001 #厚さ
    freq = 100 #周波数
    lofreq = 0 #周波数の下限
    upfreq = 0 #周波数の上限
    sampfreq = 8000.0 #サンプリング周波数
    point = 8192.0 #ポイント数
    maxmoden = 0
    mesh = 0.2
    ithq = 500
    converge = 1.e-6
    damp = zeros(0) #ダンピング
    br = zeros((0, 0))
    bz = zeros((0, 0))
    bigflux = zeros((0, 0), dtype=complex)
    bigpsz = zeros((0, 0))
    bigarea = zeros(0)
    impr = zeros(0)
    impq = zeros(0)
    impz = zeros(0)
    pr = zeros(0)
    pq = zeros(0)
    pz = zeros(0)
    area = zeros(0)
    flux = zeros((0, 0))
    conec = zeros((0, 0))
    psz = zeros((0, 0))
    #４点ガウス積分の積分点と重み係数
    zt4 = array([0.8611363116, 0.3399810436])
    wt4 = array([0.3478548451, 0.6521451549])
    
    def __init__(self):
        self.vel()  
    
    def vel(self):
        #quasi longitudinal wave velocity
        self.cb = sqrt(self.young / self.row)
        #pure longitudinal wave velocity(i.e. lateral strain is restrained)
        self.ca = sqrt(self.young * (1 - self.poisson) / (self.row * (1 + self.poisson) * (1 - 2 * self.poisson)))
        #pure longitudinal wave velocity in a plate(i.e. plain stress condition)
        self.cl = sqrt(self.young / (self.row * (1 - self.poisson * self.poisson)))
        #transversal(shear) wave velocity
        self.cs = sqrt(self.young / (self.row * 2 * (1 + self.poisson)))
        #flexural wave velocity in a beam
        self.cx = sqrt((self.thick * self.cb) * self.freq * 2 * math.pi / sqrt(12))
        #flexural wave velocity in a plate
        self.cf = sqrt((self.thick + self.cl) * self.freq * 2 * math.pi / sqrt(12))
        self.ramda = self.cx / self.freq
    
    def setparameter(self):
        while True:
            print '\n*******************************************************'
            print '\n* (P)pressure frequency --------------> %d [Hz]' % self.freq
            print '\n* (I)low and high frequency for ir --------------> from %d to %d [Hz]' % (self.lofreq, self.upfreq)
            print '\n* (S)sampling frequency --------------> %d [Hz]' % self.sampfreq
            print '\n* (L)impulse responce length --------------> %d' % self.point
            print '\n* (N)max mode number N in q direction --------------> %d' % self.maxmoden
            print '\n* (M)mesh accuracy --------------> %f' % self.mesh
            print '\n* (J)iteration max of eigen value calc --------------> %d' % self.ithq
            print '\n* (T)tolerance for checking convergence --------------> %f' % self.converge
            print '\n* (A)density of shell --------------> %f [kg/m3]' % self.row
            print '\n* (Y)Young\'s modulus --------------> %f [Pa]' % self.young
            print '\n* (R)Poisson\'s ratio --------------> ', self.poisson
            print '\n*    quasi longitudinal wave velocity'
            print '\n*       CB = %f [m/s]' % self.cb
            print '\n*    pure longitudinal wave velocity(i.e. lateral strain is restrained)'
            print '\n*       CA = %f [m/s]' % self.ca
            print '\n*    pure longitudinal wave velocity in a plate(i.e. plain stress condition)'
            print '\n*       CL = %f [m/s]' % self.cl
            print '\n*    transversal(shear) wave velocity'
            print '\n*       CS = %f [m/s]' % self.cs
            print '\n*    flexural wave velocity in a beam'
            print '\n*       CX = %f [m/s]' % self.cx
            print '\n*    flexural wave velocity in a plate'
            print '\n*       CF = %f [m/s]' % self.cf
            print '\n*    wave length for CX = %f [m]' % self.ramda
            print '\n* (H)thickness of shell --------------> %f [m]' % self.thick
            print '\n* (D)number of damping ratio --------------> ', self.damp.size
            for i in xrange(self.damp.size):
                print '\n* (%d)self.damping ratio No.%d --------------> %f' % (i+1, i+1, self.damp[i])
            print '\n*******************************************************\n'
            str = raw_input('command >>> ')
            if str[0] == 'P':
                self.freq = input('pressure frequency >>> ')
                self.vel()
            elif str[0] == 'I':
                self.lofreq = float(raw_input('low frequency for impulse responce >>> '))
                self.upfreq = float(raw_input('high frequency for impulse responce >>> '))
            elif str[0] == 'S':
                self.sampfreq = float(raw_input('sampling frequency >>> '))
            elif str[0] == 'L':
                self.point = float(raw_input('impulse response length >>> '))
            elif str[0] == 'N':
                self.maxmoden = float(raw_input('max mode number N in q direction >>> '))
            elif str[0] == 'M':
                self.mesh = float(raw_input('self.mesh accuracy >>> '))
            elif str[0] == 'J':
                self.ithq = float(raw_input('iteration max of hq >>> '))
            elif str[0] == 'T':
                self.converge = float(raw_input('tolerance for checking self.convergence >>> '))
            elif str[0] == 'A':
                self.row = float(raw_input('density of shell >>> '))
                self.vel()
            elif str[0] == 'Y':
                self.young = float(raw_input('Young\'s modulus >>> '))
                self.vel()
            elif str[0] == 'R':
                self.poisson = float(raw_input('Poisson\'s ratio >>> '))
                self.vel()
            elif str[0] == 'H':
                self.thick = float(raw_input('thickness of shell >>> '))
                self.vel()
            elif str[0] == 'D':
                tmp = float(raw_input('number of self.damping ratio >>> '))
                self.damp = zeros(tmp)
            elif str[0] == 'Q':
                break
            else:
                try:
                    i = int(str[0])
                except:
                    print 'unknown'
                    continue
                if i >= 1 and i <= self.damp.size:
                    self.damp[i - 1] = float(raw_input('self.damping ratio >>> '))
                else:
                    print 'unknown'
    
    def setshape(self):
        while True:
            print '\n*******************************************************'
            print '\n* (E)number of big element  --------------> ', size(self.br, 0)
            print '\n* (B)big element...'
            print '\n* (L)big element load file'
            print '\n*******************************************************\n'
            str = raw_input('command >>> ')
            if str[0] == 'E':
                tmp = float(raw_input('number of big element >>> '))
                self.br = zeros((tmp, 2))
                self.bz = zeros((tmp, 2))
                self.bigpsz = zeros((tmp, 2))
                self.bigflux = zeros((tmp, 3), dtype = complex)
                self.bigarea = zeros(tmp)
            elif str[0] == 'B':
                self.selectelement()
            elif str[0] == 'L':
                filename = raw_input('big element file name >>> ')
                self.loadbefile(filename)
            elif str[0] == 'Q':
                break
            else:
                print 'unknown'
    
    def selectelement(self):
        while True:
            print '\n*******************************************************'
            for i in xrange(size(self.br, 0)):
                print '\n* (%d)big element No.%d' % (i + 1, i + 1)
            print '\n*******************************************************\n'
            str = raw_input('command >>> ')
            if str[0] == 'Q':
                break
            else:
                try:
                    i = int(str)
                except:
                    print 'unknown'
                    continue
                if i >= 1 and i <=size(self.br, 0):
                    self.setelement(i - 1)
                else:
                    print 'unknown'
    
    def setelement(self, n):
        while True:
            print '\n*******************************************************'
            print '\n* (1)position 1 (r, z) --------------> (', self.br[n][0], ',', self.bz[n][0], ')'
            print '\n* (2)position 2 (r, z) --------------> (', self.br[n][1], ',', self.bz[n][1], ')'
            print '\n* (T)Type --------------> ',
            if self.bigpsz[n][0] == self.free:
                print 'Free'
            elif self.bigpsz[n][0] == self.ssupport:
                print 'Simply Supported(i.e. v=w=0), Position', self.bigpsz[n][1]
            elif self.bigpsz[n][0] == self.psupport:
                print 'Perfectly Supported(i.e. u=v=w=0), Position', self.bigpsz[n][1]
            elif self.bigpsz[n][0] == self.pointload:
                print 'Point Load, Position', self.bigpsz[n][1]
                print '\n* (Z)local z direction(from node 1 to node 2) Load -------------->', self.bigflux[n][0]
                print '\n* (S)q direction Load -------------->', self.bigflux[n][1]
                print '\n* (R)local r direction(normal direction) Load -------------->', self.bigflux[n][3]
            elif self.bigpsz[n][0] == self.lineload:
                if self.bigpsz[n][1] == 3:
                    print 'Line Load(longitudinal direction)'
                else:
                    print 'Line Load(tangential direction), Position', self.bigpsz[n][1]
                    print '\n* (A)area of line load in q direction --------------> %f ~ %f [degree]' % (-self.bigarea[n], self.bigarea[n])
                print '\n* (Z)local z direction(from node 1 to node 2) Load -------------->', self.bigflux[n][0]
                print '\n* (S)q direction Load -------------->', self.bigflux[n][1]
                print '\n* (R)local r direction(normal direction) Load -------------->', self.bigflux[n][3]
            elif self.bigpsz[n][0] == self.surfaceload:
                print 'Surface Load'
                print '\n* (A)area of surface load in q direction --------------> %f ~ %f [degree]' % (-self.bigarea[n], self.bigarea[n])
                print '\n* (Z)local z direction(from node 1 to node 2) Load -------------->', self.bigflux[n][0]
                print '\n* (S)q direction Load -------------->', self.bigflux[n][1]
                print '\n* (R)local r direction(normal direction) Load -------------->', self.bigflux[n][3]
            elif self.bigpsz[n][0] == self.axisymln:
                print 'Axisymmetric Line Load, Position', self.bigpsz[n][1]
                print '\n* (Z)local z direction(from node 1 to node 2) Load -------------->', self.bigflux[n][0]
                print '\n* (S)q direction Load -------------->', self.bigflux[n][1]
                print '\n* (R)local r direction(normal direction) Load -------------->', self.bigflux[n][2]
            elif self.bigpsz[n][0] == self.axisymsf:
                print 'Axisymmetric Surface Load'
                print '\n* (Z)local z direction(from node 1 to node 2) Load -------------->', self.bigflux[n][0]
                print '\n* (S)q direction Load -------------->', self.bigflux[n][1]
                print '\n* (R)local r direction(normal direction) Load -------------->', self.bigflux[n][3]
            print '\n*******************************************************\n'
            str = raw_input('command >>> ')
            if str[0] == '1':
                self.br[n][0] = float(raw_input('position1 r >>> '))
                self.bz[n][0] = float(raw_input('position1 z >>> '))
            elif str[0] == '2':
                self.br[n][1] = float(raw_input('position2 r >>> '))
                self.bz[n][1] = float(raw_input('position2 z >>> '))
            elif str[0] == 'T':
                self.bigpsz[n][0] = float(raw_input('Type(0:Free, 1:Simply Supported, 3:Point Load, 4:Line Load, 5:Surface Load,\n     6:Axisymmetric Line Load, 7:Axisymmetric Surface Load, 8:Perfectly Supported) >>> '))
                if self.bigpsz[n][0] != self.free and self.bigpsz[n][0] != self.surfaceload and self.bigpsz[n][0] != self.axisymsf:
                    self.bigpsz[n][1] = float(raw_input('Position >>> '))
            elif str[0] == 'A':
                self.bigarea[n] = float(raw_input('area of load in q direction[degree] >>> '))
            elif str[0] == 'Z':
                tmpreal = float(raw_input('local z Load real >>> '))
                tmpimag = float(raw_input('local z Load imag >>> '))
                self.bigflux[n][0] = complex(tmpreal, tmpimag)
            elif str[0] == 'S':
                tmpreal = float(raw_input('q Load real >>> '))
                tmpimag = float(raw_input('q Load imag >>> '))
                self.bigflux[n][1] = complex(tmpreal, tmpimag)
            elif str[0] == 'R':
                tmpreal = float(raw_input('local r Load real >>> '))
                tmpimag = float(raw_input('local r Load imag >>> '))
                self.bigflux[n][2] = complex(tmpreal, tmpimag)
            elif str[0] == 'Q':
                break
            else:
                print 'unknown'
    
    def setlisten(self):
        while True:
            print '\n*******************************************************'
            print '\n* (I)number of impulse response listening point -------------->', self.impr.size
            print '\n* (M)impulse response listening point ...'
            print '\n* (F)impulse response listening point load file'
            print '\n* (P)number of displacement listening point  -------------->', self.pr.size
            print '\n* (R)displacement listening point ...'
            print '\n* (E)displacement listening point load file'
            print '\n*******************************************************\n'
            str = raw_input('command >>> ')
            if str[0] == 'I':
                nu = float(raw_input('number of impulse response listening point >>> '))
                if nu > 10:
                    nu = 10
                self.impr = zeros(nu)
                self.impq = zeros(nu)
                self.impz = zeros(nu)
            elif str[0] == 'M':
                self.setimptau()
            elif str[0] == 'F':
                filename = raw_input('impulse response listening point load file name >>> ')
                self.loadlpfile(filename)
            elif str[0] == 'P':
                nu = float(raw_input('number of displacement listening point >>> '))
                self.pr = zeros(nu)
                self.pq = zeros(nu)
                self.pz = zeros(nu)
            elif str[0] == 'R':
                self.setprestau()
            elif str[0] == 'E':
                filename = raw_input('displacement listening point load file >>> ')
                self.loadlpfile(filename)
            elif str[0] == 'Q':
                break
            else:
                print 'unknown'
    
    def setimptau(self):
        while True:
            print '\n*******************************************************'
            if self.impr.size > 0:
                print '\n* (1)impulse response listening point No.1 (r[m], q[degree], z[m]) --------------> (', self.impr[0], ',', self.impq[0], ',', self.impz[0], ')'
            if self.impr.size > 1:
                print '\n*   :                :                  :                  :                 :'
                print '\n* (', self.impr.size, 'impulse response listening self.point No.', self.impr.size, '(r[m], q[degree], z[m]) --------------> (', self.impr[self.impr.size - 1], ',', self.impq[self.impr.size - 1], ',', self.impz[self.impr.size - 1], ')'
            print '\n*******************************************************\n'
            str = raw_input('command >>> ')
            if str[0] == 'Q':
                for i in xrange(self.impr.size):
                    print '(', i, ')impulse response listening point No.', i, '(r[m], q[degree], z[m]) --------------> (', self.impr[i - 1], ',', self.impq[i - 1], ',', self.impz[i - 1], ')'
                break
            else:
                try:
                    i = int(str)
                except:
                    print 'unknown'
                    continue
                if i >= 1 and i <= self.impr.size:
                    self.impr[i - 1] = float(raw_input('listening point r[m] >>> '))
                    self.impq[i - 1] = float(raw_input('listening point q[degree] >>> '))
                    self.impz[i - 1] = float(raw_input('listening point z[m] >>> '))
                    print '(', i, ')impulse response listening point No.', i, '(r[m], q[degree], z[m]) --------------> (', self.impr[i - 1], ',', self.impq[i - 1], ',', self.impz[i - 1], ')'
                else:
                    print 'unknown'
    
    
    def setprestau(self):
        while True:
            print '\n*******************************************************'
            if self.pr.size > 0:
                print '\n* (1)displacement listening point No.1 (r[m], q[degree], z[m]) --------------> (', self.pr[0], ',', self.pq[0], ',', self.pz[0], ')'
            if self.pr.size > 1:
                print '\n*   :                :                  :                  :                 :'
                print '\n* (', self.pr.size, 'displacement listening point No.', self.pr.size, '(r[m], q[degree], z[m]) --------------> (', self.pr[self.pr.size - 1], ',', self.pq[self.pr.size - 1], ',', self.pz[self.pr.size - 1], ')'
            print '\n*******************************************************\n'
            str = raw_input('command >>> ')
            if str[0] == 'Q':
                for i in xrange(self.pr.size):
                    print '(', i, ')displacement listening point No.', i, '(r[m], q[degree], z[m]) --------------> (', self.pr[i - 1], ',', self.pq[i - 1], ',', self.pz[i - 1], ')'
                break
            else:
                try:
                    i = int(str)
                except:
                    print 'unknown'
                    continue
                if i >= 1 and i <= self.pr.size:
                    self.pr[i - 1] = float(raw_input('listening point r[m] >>> '))
                    self.pq[i - 1] = float(raw_input('listening point q[degree] >>> '))
                    self.pz[i - 1] = float(raw_input('listening point z[m] >>> '))
                    print '(', i, ')displacement listening point No.', i, '(r[m], q[degree], z[m]) --------------> (', self.pr[i - 1], ',', self.pq[i - 1], ',', self.pz[i - 1], ')'
                else:
                    print 'unknown'
    
    def calc(self):
        while True:
            print '\n*******************************************************'
            print '\n* (N)calculate natural frequencies'
            print '\n* (F)print natural frequencies'
            print '\n* (R)write to file element condition and node coordinates'
            print '\n* (P)displacement distribution...'
            print '\n* (I)impulse response...'
            print '\n*******************************************************\n'
            str = raw_input('command >>> ')
            if str[0] == 'N':
                cmode = float(raw_input('circumferential mode number'))
                print '########## calculate natural frequencies ##########'
                print '########## check divide number ##########'
                divide, eps = self.checkdiv()
                print ' divide =', divide, ' eps =', eps
                print '########## cut self.mesh ##########'
                noder, nodez = self.dgnode(divide, eps)
                self.mkelement(divide, eps, noder, nodez, Hsp, speaker)
                print ' node =', noder.size
                print ' number of element =', size(self.conec, 0)
                if cmode == 0:
                    print '########## make M and K matrix for mode 0 ##########'
                    self.makeMCKf0(MM, CC, KK, FF, evalue, evector, self.damp, noder, nodez, self.conec, psz, flux, area,
                                  True, False, True, True, self.ithq)
                    print '########## apply boundary condition for mode 0 ##########'
                    rd = self.reduceBV0(MM, KK, FF, self.conec, psz)
                else:
                    print '########## make M and K matrix for mode', cmode, '##########'
                    self.makeMCKf(MM, CC, KK, FF, evalue, evector, self.damp, noder, nodez, self.conec, psz, flux, area, cmode,
                                 True, False, True, True, self.ithq)
                    print '########## apply boundary condition for mode', cmode, '##########'
                    rd = self.reduceBV(MM, KK, FF, self.conec, psz)
                print '########## compute eigenvalue by hq ##########'
                evalue = zeros(size(KK, 0) - rd)
                evector = zeros((size(KK, 0), size(KK, 0)))
                gshqlB(KK, MM, KK.row()-rd, evalue, evector, self.ithq)
            elif str[0] == 'F':
                pmode = float(raw_input('print max mode number >>> '))
                if pmode > evalue.size:
                    pmode = evalue.size
                print 'mode number\tnatural frequencies[Hz]'
                for i in xrange(pmode):
                    f = sqrt(evalue[i]) / 2 * math.pi
                    print i, '\t', f
            elif str[0] == 'R':
                print '########## check divide number ##########'
                divide, eps = self.checkdiv()
                print ' divide =', divide, ' eps =', eps
                print '########## cut self.mesh ##########'
                noder, nodez = self.dgnode(divide, eps)
                self.mkelement(divide, eps, noder, nodez, Hsp, speaker)
                print ' node =', noder.size
                print ' number of element =', size(self.conec, 0)
                filename = raw_input('element condition file name >>> ')
                self.elementsave(filename, noder, nodez)
                filename = raw_input('node coordinates file name >>> ')
                self.nodesave(filename, noder, nodez)
            elif str[0] == 'P':
                self.pressure()
            elif str[0] == 'I':
                self.impulse()
            elif str[0] == 'Q':
                break
            else:
                print 'unknown'
    
    def pressure(self):
        filename = ''
        while True:
            print '\n*******************************************************'
            print '\n* (F)file name -------------->', filename
            print '\n* (S)start'
            print '\n*******************************************************\n'
            str = raw_input('command >>> ')
            if str[0] == 'F':
                filename = raw_input('file name >>> ')
            elif str[0] == 'S':
                uu = zeros((1, self.pr.size), dtype = complex)
                vv = zeros((1, self.pr.size), dtype = complex)
                ww = zeros((1, self.pr.size), dtype = complex)
                bb = zeros((1, self.pr.size), dtype = complex)
                print '########## format data save ##########'
                self.dmyuvwbsave(filename, self.pr.size)
                omega = freq * 2 * math.pi
                print '########## check divide number ##########'
                divide, eps = self.checkdiv()
                print ' divide =', divide, ' eps =', eps
                print '########## cut self.mesh ##########'
                noder, nodez = self.dgnode(divide, eps)
                self.mkelement(divide, eps, noder, nodez, Hsp, speaker)
                print ' node =', noder.size
                print ' number of element =', size(self.conec, 0)
                print '########## make matrix for mode 0 ##########'
                B, V = self.makeBV0(noder, nodez, omega)
                print '########## apply boundary condition for mode 0 ##########'
                self.applybc0(B, V)
                print '########## solve matrix by srfcho ##########'
                srfcho(B, size(B, 0))
                srfchosv(B, PP, size(B, 0))
                print '########## calculate displacement for mode 0 ##########'
                for i in xrange(self.pr.size):
                    un, wn, bn = self.calcuwb0(self.pr[i], self.pz[i], noder, nodez, V, eps)
                    uu[i] = un
                    vv[i] = 0
                    ww[i] = wn
                    bb[i] = bn
                    print ' listening point (', self.pr[i], ',', self.pq[i], ',', self.pz[i], ')'
                    print '                 un =', un
                    print '                 wn =', wn
                    print '                 bn =', bn
                    self.wrtfile(filename, i, uu[i], vv[i], ww[i], bb[i])
                for i in xrange(self.maxmoden):
                    print '########## make matrix for mode', i + 1, '##########'
                    self.makeBV(B, PP, noder, nodez, self.conec, psz, flux, area, i + 1, omega)
                    print '########## apply boundary condition for mode', i + 1, '##########'
                    self.applybc(B, V, self.conec, psz)
                    print '########## solve matrix by srfcho ##########'
                    srfcho(B, size(B, 0))
                    srfchosv(B, V, size(B, 0))
                    print '########## calculate displacement for mode', i + 1, '##########'
                    for j in xrange(pr.size):
                        self.calcuvwb(un, vn, wn, bn, pr[j], pz[j], noder, nodez, self.conec, V, eps)
                        uu[j] += un * cos((i + 1) * pq[j] * math.pi / 180)
                        vv[j] += vn * cos((i + 1) * pq[j] * math.pi / 180)
                        ww[j] += wn * cos((i + 1) * pq[j] * math.pi / 180)
                        bb[j] += bn * cos((i + 1) * pq[j] * math.pi / 180)
                        print ' listening point (', pr[j], ',', pq[j], ',', pz[j], ')'
                        print '                 un =', un
                        print '                 vn =', vn
                        print '                 wn =', wn
                        print '                 bn =', bn
                        fil.wrtfile(filename, j, uu[j], vv[j], ww[j], bb[j])
            elif str[0] == 'Q':
                break
            else:
                print 'unknown'
    
    def impulse(self):
        filename = ''
        remesh, speaker = True, False
        A = 5.62
        Ld = 0.97e-3
        Red = 7.4
        while True:
            print '\n*******************************************************'
            print '\n* (F)file name -------------->', filename
            print '\n* (R)remesh -------------->',
            if remesh:
                print 'on'
            else:
                print 'off'
            print '\n* (D)speaker --------------->',
            if speaker:
                print 'on'
            else:
                print 'off'
            print '\n* (S)start'
            if speaker:
                print '\n (A)力係数 --------------->', A
                print '\n (L)コイルのインダクタンス --------------->', Ld
                print '\n (RE)電気抵抗 --------------->', Red
            print '\n*******************************************************\n'
            str = raw_input('command >>> ')
            if str[0] == 'F':
                filename = raw_input('file name >>> ')
                print 'displacement is saved as %s_*.png' % filename
            elif str[0] == 'R':
                remesh = float(raw_input('remesh(0:off, 1:on) >>> '))
            elif str[0] == 'D':
                speaker = float(raw_input('speaker(0:off, 1:on) >>> '))
            elif str[0] == 'A':
                A = float(raw_input(u'力係数 >>> '))
            elif str[0] == 'L':
                Ld = float(raw_input(u'コイルのインダクタンス >>> '))
            elif str[0] == 'RE':
                Red = float(raw_input(u'電気抵抗 >>> '))
            elif str[0] == 'S':
                uu = array([], dtype = complex)
                vv = array([], dtype = complex)
                ww = array([], dtype = complex)
                bb = array([], dtype = complex)
                plotn = 0
                plotfreq = array([])
                for i in xrange(self.point / 2 + 1):
                    self.freq = self.sampfreq / self.point * i
                    if self.freq >= self.lofreq and self.freq <=self.upfreq:
                        plotfreq = append(plotfreq, self.freq)
                        print '############## calculate frequency :', self.freq, '[Hz] ##############'
                        omega = self.freq * 2 * math.pi
                        if remesh:
                            self.vel()
                        if remesh or not meshdone:
                            print '########## check divide number ##########'
                            divide, eps = self.checkdiv()
                            print ' divide =', divide, ' eps =', eps
                            print '########## cut self.mesh ##########'
                            noder, nodez = self.dgnode(divide, eps)
                            jomega = complex(0, omega)
                            Hsp = A / (jomega * Ld + Red)
                            self.mkelement(divide, eps, noder, nodez, Hsp, speaker)
                            print ' node =', noder.size
                            print ' number of element =', size(self.conec, 0)
                            meshdone = True
                        print '########## make matrix for mode 0 ##########'
                        B, V = self.makeBV0(noder, nodez, omega, speaker, A)
                        print '########## apply boundary condition for mode 0 ##########'
                        self.applybc0(B, V)
                        print '########## solve matrix ##########'
                        V = solve(B, V)
                        print '########## calculate displacement for mode 0 ##########'
                        for j in xrange(self.impr.size):
                            un, wn, bn = self.calcuwb0(self.impr[j], self.impz[j], noder, nodez, V, eps)
                            uu = append(uu, un)
                            vv = append(vv, 0)
                            ww = append(ww, wn)
                            bb = append(bb, bn)
                            print ' listening point (', self.impr[j], ',', self.impq[j], ',', self.impz[j], ')'
                            print '                 un =', un
                            print '                 wn =', wn
                            print '                 bn =', bn
                        for j in xrange(self.maxmoden):
                            print '########## make matrix for mode', j + 1, '##########'
                            self.makeBV(B, V, noder, nodez, j + 1, omega)
                            print '########## apply boundary condition for mode', j + 1, '##########'
                            self.applybc(B, V)
                            print '########## solve matrix by srfcho ##########'
                            srrcho(B, size(B, 0))
                            srfchosv(B, V, size(B, 0))
                            print '########## calculate displacement for mode', j + 1, '##########'
                            for k in xrange(self.impr.size):
                                self.calcuvwb(un, vn, wn, bn, self.impr[k], self.impz[k], noder, nodez, V, eps)
                                uu[self.impr.size * plotn + k] += un * cos((j + 1) * self.impq[k] * math.pi / 180)
                                vv[self.impr.size * plotn + k] += vn * sin((j + 1) * self.impq[k] * math.pi / 180)
                                ww[self.impr.size * plotn + k] += wn * cos((j + 1) * self.impq[k] * math.pi / 180)
                                bb[self.impr.size * plotn + k] += bn * cos((j + 1) * self.impq[k] * math.pi / 180)
                                print ' listening point (', self.impr[k], ',', self.impq[k], ',', self.impz[k], ')'
                                print '                 un =', un
                                print '                 vn =', vn
                                print '                 wn =', wn
                                print '                 bn =', bn
                        plotn += 1
                self.plot(uu, vv, ww, bb, plotn, plotfreq, filename)
            elif str[0] == 'Q':
                break
            else:
                print 'unknown'
                
                
    
    #以下、入出力関係
    def loadbefile(self, filename):
        try:
            ft = file(filename)
        except:
            print 'File not found or illegal access'
            sys.exit()
        ss = 0
        while True:
            line = ft.readline()
            if line == '':
                break
            elif line[0] != '\n':
                ss += 1
        ft.seek(0)
        bne = ss / 3
        self.br = zeros((bne, 2))
        self.bz = zeros((bne, 2))
        self.bigpsz = zeros((bne, 2))
        self.bigflux = zeros((bne, 3))
        self.bigarea = zeros(bne)
        for i in xrange(bne):
            line = ft.readline()
            token = line.split()
            if token[0][0] == 'G':
                self.bigpsz[i - 1][0] = self.gsupport
                self.bigpsz[i - 1][1] = int(token[1])
                self.bigflux[i - 1][0] = complex(float(token[2]), float(token[3]))
                self.bigflux[i - 1][1] = complex(float(token[4]), float(token[5]))
                self.bigflux[i - 1][2] = complex(float(token[6]), float(token[7]))
                self.br[i - 1][0] = float(token[8])
                self.br[i - 1][1] = float(token[10])
                self.bz[i - 1][0] = float(token[9])
                self.bz[i - 1][1] = float(token[11])
            elif token[0][0] == 'S':
                self.bigpsz[i - 1][0] = self.ssupport
                self.bigpsz[i - 1][1] = int(token[1])
                self.br[i - 1][0] = float(token[2])
                self.br[i - 1][1] = float(token[4])
                self.bz[i - 1][0] = float(token[3])
                self.bz[i - 1][1] = float(token[5])
            elif token[0][0] == 'T':
                self.bigpsz[i - 1][0] = self.psupport
                self.bigpsz[i - 1][1] = int(token[1])
                self.br[i - 1][0] = float(token[2])
                self.br[i - 1][1] = float(token[4])
                self.bz[i - 1][0] = float(token[3])
                self.bz[i - 1][1] = float(token[5])
            elif token[0][0] == 'F':
                self.bigpsz[i - 1][0] = self.free
                self.br[i - 1][0] = float(token[1])
                self.br[i - 1][1] = float(token[3])
                self.bz[i - 1][0] = float(token[2])
                self.bz[i - 1][1] = float(token[4])
            elif token[0][0] == 'P':
                self.bigpsz[i - 1][0] = self.pointload
                self.bigpsz[i - 1][1] = int(token[1])
                self.bigflux[i - 1][0] = complex(float(token[2]), float(token[3]))
                self.bigflux[i - 1][1] = complex(float(token[4]), float(token[5]))
                self.bigflux[i - 1][2] = complex(float(token[6]), float(token[7]))
                self.br[i - 1][0] = float(token[8])
                self.br[i - 1][1] = float(token[10])
                self.bz[i - 1][0] = float(token[9])
                self.bz[i - 1][1] = float(token[11])
            elif token[0][0] == 'L':
                self.bigpsz[i - 1][0] = self.lineload
                self.bigpsz[i - 1][1] = int(token[1])
                self.bigarea[i - 1] = float(token[2])
                self.bigflux[i - 1][0] = complex(float(token[3]), float(token[4]))
                self.bigflux[i - 1][1] = complex(float(token[5]), float(token[6]))
                self.bigflux[i - 1][2] = complex(float(token[7]), float(token[8]))
                self.br[i - 1][0] = float(token[9])
                self.br[i - 1][1] = float(token[11])
                self.bz[i - 1][0] = float(token[10])
                self.bz[i - 1][1] = float(token[12])
            elif token[0][0] == 'A':
                self.bigpsz[i - 1][0] = self.surfaceload
                self.bigarea[i - 1] = float(token[1])
                self.bigflux[i - 1][0] = complex(float(token[2]), float(token[3]))
                self.bigflux[i - 1][1] = complex(float(token[4]), float(token[5]))
                self.bigflux[i - 1][2] = complex(float(token[6]), float(token[7]))
                self.br[i - 1][0] = float(token[8])
                self.br[i - 1][1] = float(token[10])
                self.bz[i - 1][0] = float(token[9])
                self.bz[i - 1][1] = float(token[11])
            elif token[0][0] == 'Y':
                self.bigpsz[i - 1][0] = self.axisymln
                self.bigpsz[i - 1][1] = int(token[1])
                self.bigflux[i - 1][0] = complex(float(token[2]), float(token[3]))
                self.bigflux[i - 1][1] = complex(float(token[4]), float(token[5]))
                self.bigflux[i - 1][2] = complex(float(token[6]), float(token[7]))
                self.br[i - 1][0] = float(token[8])
                self.br[i - 1][1] = float(token[10])
                self.bz[i - 1][0] = float(token[9])
                self.bz[i - 1][1] = float(token[11])
            elif token[0][0] == 'M':
                self.bigpsz[i - 1][0] = self.axisymsf
                self.bigflux[i - 1][0] = complex(float(token[1]), float(token[2]))
                self.bigflux[i - 1][1] = complex(float(token[3]), float(token[4]))
                self.bigflux[i - 1][2] = complex(float(token[5]), float(token[6]))
                self.br[i - 1][0] = float(token[7])
                self.br[i - 1][1] = float(token[9])
                self.bz[i - 1][0] = float(token[8])
                self.bz[i - 1][1] = float(token[10])
            else:
                print 'big element file error'
                sys.exit()
        ft.close()
    
    def loadlpfile(self, filename):
        try:
            ft = file(filename)
        except:
            print 'File not found or illegal access'
            sys.exit()
        ss = 0
        while True:
            line = ft.readline()
            if line == '':
                break
            elif line[0] != '\n':
                ss += 1
        ft.seek(0)
        self.impr = zeros(ss)
        self.impq = zeros(ss)
        self.impz = zeros(ss)
        for i in xrange(ss):
            line = ft.readline()
            token = line.split()
            self.impr[i - 1] = float(token[0])
            self.impq[i - 1] = float(token[1])
            self.impz[i - 1] = float(token[2])
        ft.close()
    
    def elementsave(self, filename, noder, nodez):
        pass
    
    def nodesave(self, filename, noder, nodez):
        pass
    
    def dmyuvwbsave(self, filename, sleng):
        pass
    
    def wrtfile(self, filename, j, wu, wv, ww, wb):
        pass
    
    def wrtirfile(self, tmpname, i, wu):
        pass
    
    def plot(self, uu, vv, ww, bb, plotn, plotfreq, filename):
        for i in xrange(self.impr.size):
                    u = array([uu[j] for j in xrange(self.impr.size * plotn) if j % self.impr.size == i])
                    v = array([vv[j] for j in xrange(self.impr.size * plotn) if j % self.impr.size == i])
                    w = array([ww[j] for j in xrange(self.impr.size * plotn) if j % self.impr.size == i])
                    b = array([bb[j] for j in xrange(self.impr.size * plotn) if j % self.impr.size == i])
                    ##### 変位をプロット #####
                    urange = max(abs(u))
                    vrange = max(abs(v))
                    wrange = max(abs(w))
                    brange = max(abs(b))
                    pylab.clf()
                    pylab.subplot(221)
                    pylab.semilogx(plotfreq, 20 * log10(abs(u) / urange))
                    pylab.ylabel('u[dB]')
                    pylab.xlabel('frequency[Hz]')
                    pylab.grid(True)
                    pylab.subplot(222)
                    if vrange != 0:
                        pylab.semilogx(plotfreq, 20 * log10(abs(v) / vrange))
                    else:
                        pylab.semilogx(plotfreq, abs(v))
                    pylab.ylabel('v[dB]')
                    pylab.xlabel('frequency[Hz]')
                    pylab.grid(True)
                    pylab.subplot(223)
                    pylab.semilogx(plotfreq, 20 * log10(abs(w) / wrange))
                    pylab.ylabel('w[dB]')
                    pylab.xlabel('frequency[Hz]')
                    pylab.grid(True)
                    pylab.subplot(224)
                    pylab.semilogx(plotfreq, 20 * log10(abs(b) / brange))
                    pylab.ylabel('b[dB]')
                    pylab.xlabel('frequency[Hz]')
                    pylab.grid(True)
                    pylab.subplots_adjust(wspace = 0.3, hspace = 0.3)
                    pylab.savefig(filename + 'dis_%d.png' % (i + 1))
                    ##### 速度をプロット #####
                    urange = max(abs(2 * math.pi * plotfreq * u))
                    vrange = max(abs(2 * math.pi * plotfreq * v))
                    wrange = max(abs(2 * math.pi * plotfreq * w))
                    brange = max(abs(2 * math.pi * plotfreq * b))
                    pylab.clf()
                    pylab.subplot(221)
                    pylab.semilogx(plotfreq, 20 * log10(abs(2 * math.pi * plotfreq * u) / urange))
                    pylab.ylabel('u[dB]')
                    pylab.xlabel('frequency[Hz]')
                    pylab.grid(True)
                    pylab.subplot(222)
                    if vrange != 0:
                        pylab.semilogx(plotfreq, 20 * log10(abs(2 * math.pi * plotfreq * v) / vrange))
                    else:
                        pylab.semilogx(plotfreq, abs(plotfreq * v))
                    pylab.ylabel('v[dB]')
                    pylab.xlabel('frequency[Hz]')
                    pylab.grid(True)
                    pylab.subplot(223)
                    pylab.semilogx(plotfreq, 20 * log10(abs(2 * math.pi * plotfreq * w) / wrange))
                    pylab.ylabel('w[dB]')
                    pylab.xlabel('frequency[Hz]')
                    pylab.grid(True)
                    pylab.subplot(224)
                    pylab.semilogx(plotfreq, 20 * log10(abs(2 * math.pi * plotfreq * b) / brange))
                    pylab.ylabel('b[dB]')
                    pylab.xlabel('frequency[Hz]')
                    pylab.grid(True)
                    pylab.subplots_adjust(wspace = 0.3, hspace = 0.3)
                    pylab.savefig(filename + 'vel_%d.png' % (i + 1))

    
    
    #以下、ソルバー
    def checkdiv(self):
        nbe = size(self.br, 0) #big elementの数
        max = 0
        for i in xrange(nbe):
            sr = self.br[i][1] - self.br[i][0] #r2-r1
            sz = self.bz[i][1] - self.bz[i][0] #z2-z1
            sss = sr**2 + sz**2 #距離の２乗
            if i == 0:
                min = sss
            elif min > sss:
                min = sss
            if max < sss:
                max = sss
        hen = max = sqrt(max) #距離
        min = sqrt(min)
        me_ramda = self.mesh * self.ramda #波長にmesh accuracyをかけて、それをひとつの要素の長さとする
        divide = 1
        while me_ramda < hen: #big elementの変の長さを分割するのを繰り返し、me_ramdaよりも小さくする。
            divide += 1
            hen = max / divide
        eps = min / divide / 100
        print 'min. mesh =', min / divide #最小の要素長
        print 'max. mesh =', max / divide #最大の要素長
        return divide, eps
    
    def dgnode(self, divide, eps): #全節点の座標(r,z)を計算
        nbe = size(self.br, 0) #big elementの数
        div1 = divide + 1 #１つのbig elementにおける節点の数
        node = div1 * nbe #全big elementにおける節点の数(総節点数)
        ttr = zeros(node)
        ttz = zeros(node)
        gnode = 0
        bignode = 0
        for i in xrange(nbe):
            b0r = self.br[i][0]
            b0z = self.bz[i][0]
            sr = (self.br[i][1] - b0r) / divide
            sz = (self.bz[i][1] - b0z) / divide
            for j in xrange(div1):
                tmpr = b0r + sr * j
                tmpz = b0z + sz * j
                cores = False
                for k in xrange(bignode):
                    if abs(ttr[k] - tmpr) <= eps and abs(ttz[k] - tmpz) <= eps:
                        cores = True
                        break
                if cores == False:
                    gnode += 1
                    ttr[gnode - 1] = tmpr #節点のr座標
                    ttz[gnode - 1] = tmpz #節点のz座標
            bignode = gnode
        noder = zeros(gnode)
        nodez = zeros(gnode)
        for i in xrange(gnode):
            noder[i] = ttr[i]
            nodez[i] = ttz[i]
        return noder, nodez
    
    def mkelement(self, divide, eps, noder, nodez, Hsp, speaker):
        nbe = size(self.br, 0) #big elementの数
        node = noder.size #節点の数
        ne = divide * nbe #要素の数
        self.conec = zeros((ne, 2))
        self.psz = zeros((ne, 2))
        self.flux = zeros((ne, 3))
        self.area = zeros(ne)
        nne = 0
        for i in xrange(nbe): #loop for big element
            type = self.bigpsz[i][0] #big elementの境界条件のtype
            where = self.bigpsz[i][1] #big elementの境界条件の位置
            b0r = self.br[i][0] #big elementのr1座標
            b0z = self.bz[i][0] #big elementのz1座標
            sr = (self.br[i][1] - b0r) / divide #要素一つのr方向の長さ
            sz = (self.bz[i][1] - b0z) / divide #要素一つのz方向の長さ
            for j in xrange(divide): #make element per one big element
                nne += 1
                p1r = b0r + sr * j #j+1番目の要素の一つ目の節点のr座標
                p1z = b0z + sz * j #j+1番目の要素の一つ目の節点のz座標
                p2r = b0r + sr * (j + 1) #j+1番目の要素の二つ目の節点のr座標
                p2z = b0z + sz * (j + 1) #j+1番目の要素の二つ目の節点のz座標
                cores1 = False
                cores2 = False
                for k in xrange(node): #節点に番号を振っていく[1,2;3,4;5,6;…]
                    if cores1 == False:
                        if abs(noder[k] - p1r) <= eps and abs(nodez[k] - p1z) <= eps:
                            self.conec[nne - 1][0] = k + 1
                            cores1 = True
                    if cores2 == False:
                        if abs(noder[k] - p2r) <= eps and abs(nodez[k] - p2z) <= eps:
                            self.conec[nne - 1][1] = k + 1
                            cores2 = True
                    if cores1 == True and cores2 == True:
                        break
                if cores1 == False or cores2 == False:
                    print 'making element error'
                    sys.exit()
                if type == self.ssupport or type == self.psupport or type == self.lineload or type == self.axisymln:
                    if where == 1:
                        if j == 0:
                            self.psz[nne - 1][0] = type
                            self.psz[nne - 1][1] = 1
                            if speaker:
                                self.flux[nne - 1] = self.bigflux[i] * Hsp
                            else:
                                self.flux[nne - 1] = self.bigflux[i]
                            self.area[nne - 1] = self.bigarea[i]
                        else:
                            self.psz[nne - 1][0] = self.free
                    elif where == 2:
                        if j == divide - 1:
                            self.psz[nne - 1][0] = type
                            self.psz[nne - 1][1] = 2
                            if speaker:
                                self.flux[nne - 1] = self.bigflux[i] * Hsp
                            else:
                                self.flux[nne - 1] = self.bigflux[i]
                            self.area[nne - 1] = self.bigarea[i]
                        else:
                            self.psz[nne - 1][0] = self.free
                    elif where == 3:
                        self.psz[nne - 1][0] = type
                        self.psz[nne - 1][1] = 3
                        self.flux[nne - 1] =self.bigflux[i]
                    else:
                        print 'Simply Supported or Perfectly Supported or (Axisymmetric) Line load set error'
                        sys.exit()
                elif type == self.pointload:
                    if where == 1:
                        if j == 0:
                            self.psz[nne - 1][0] = type
                            self.psz[nne - 1][1] = 1
                            self.flux[nne - 1] = self.bigflux[i]
                        else:
                            self.psz[nne - 1][0] = self.free
                    elif where == 2:
                        if j == divide - 1:
                            self.psz[nne - 1][0] = type
                            self.psz[nne - 1][1] = 2
                            self.flux[nne - 1] = self.bigflux[i]
                        else:
                            self.psz[nne - 1][0] = self.free
                    else:
                        print 'Point Load set error'
                        sys.exit()
                elif type == self.free or type == self.surfaceload or type == self.axisymsf:
                    self.psz[nne - 1][0] = type
                    self.flux[nne - 1] = self.bigflux[i]
                    self.area[nne - 1] = self.bigarea[i]
                else:
                    print 'Big Element Type set error'
                    sys.exit()
    
    def makeMCKf0(self, MM, CC, KK, FF, evalue, evector, noder, nodez, calcMM, caclCC, calcKK, calcff):
        ne = size(self.conec, 0)
        mnode = noder.size
        mnode3 = mnode * 3
        Nn1 = zeros((2, 3))
        Nn2 = zeros((2, 3))
        Bn1 = zeros((4, 3))
        Bn2 = zeros((4, 3))
        DD = zeros((4, 4))
        tmp1 = zeros((3, 3))
        smp1 = zeros((3, 3))
        ans11 = zeros((3, 3))
        ans12 = zeros((3, 3))
        ans22 = zeros((3, 3))
        ansff = zeros((3, 2))
        self.ramda = zeros((3, 3))
    
    
    def reduceBV0(self, MM, KK, FF):
        pass
    
    def makeMCKf(self, MM, CC, KK, FF, evalue, evector, noder, nodez, cmode, calcMM, calcCC, calcKK, calcff):
        pass
    
    def reduceBV(self, MM, KK, FF):
        pass
    
    def makeBV0(self, noder, nodez, omega, speaker, A):
        jomega = complex(0, omega)
        ne = size(self.conec, 0) #要素数
        mnode = noder.size #節点数
        mnode3 = mnode * 3
        B = zeros((mnode3, mnode3), dtype = complex)
        V = zeros(mnode3, dtype = complex)
        OMOM = omega**2
        N01 = zeros((2, 3))
        N02 = zeros((2, 3))
        B01 = zeros((4, 3))
        B02 = zeros((4, 3))
        DD = zeros((4, 4)) #弾性マトリクス作成
        tmp1 = zeros((3, 3))
        ans11 = zeros((3, 3))
        ans12 = zeros((3, 3))
        ans22 = zeros((3, 3))
        bns11 = zeros((3, 3))
        bns12 = zeros((3, 3))
        bns22 = zeros((3, 3))
        ramda = zeros((3, 3))
        ##############################弾性マトリクス作成##############################
        DD[0][0] = DD[1][1] = self.young * self.thick / (1 - self.poisson**2)
        DD[1][0] = DD[0][1] = DD[0][0] * self.poisson
        DD[2][2] = DD[3][3] = self.thick**2 * DD[0][0] / 12
        DD[3][2] = DD[2][3] = self.thick**2 * DD[0][1] / 12
        #########################################################################
        for i in xrange(ne): #loop for elements
            g1 = self.conec[i][0]
            g2 = self.conec[i][1]
            st1 = (g1 - 1) * 3
            st2 = (g2 - 1) * 3
            u1 = noder[g2 - 1] - noder[g1 - 1] #要素のr方向の長さ
            v1 = nodez[g2 - 1] - nodez[g1 - 1] #要素のz方向の長さ
            L = sqrt(u1**2 + v1**2)
            cf = v1 / L #cosΦ
            sf = u1 / L #sinΦ
            ##########マトリクスλ作成##########
            ramda[0][0] = cf
            ramda[0][1] = sf
            ramda[1][0] = -sf
            ramda[1][1] = cf
            ramda[2][2] = 1
            ###############################
            integ11 = zeros((3, 3))
            integ12 = zeros((3, 3))
            integ22 = zeros((3, 3))
            hnteg11 = zeros((3, 3))
            hnteg12 = zeros((3, 3))
            hnteg22 = zeros((3, 3))
            for j in xrange(2): #numerical integral loop
                #一つ目のxiにおいて計算
                xi = (1 - self.zt4[j]) / 2 #区間(0,1)で４点ガウス積分する際の積分点
                rr = noder[g1 - 1] + (noder[g2 - 1] - noder[g1 - 1]) * xi #数値積分する際のr座標
                self.setN0B0(N01, N02, B01, B02, xi, L, rr, cf, sf) #形状マトリクスNとひずみマトリクスBを計算
                tmp1 = rr * self.wt4[j] * dot(transpose(N01), N01)
                integ11 += tmp1
                tmp1 = rr * self.wt4[j] * dot(transpose(N01), N02)
                integ12 += tmp1
                tmp1 = rr * self.wt4[j] * dot(transpose(N02), N02)
                integ22 += tmp1
                tmp1 = rr * self.wt4[j] * dot(dot(transpose(B01), DD), B01)
                hnteg11 += tmp1
                tmp1 = rr * self.wt4[j] * dot(dot(transpose(B01), DD), B02)
                hnteg12 += tmp1
                tmp1 = rr * self.wt4[j] * dot(dot(transpose(B02), DD), B02)
                hnteg22 += tmp1
                #二つ目のxiにおいて計算
                xi = (1 + self.zt4[j]) / 2
                rr = noder[g1 - 1] + (noder[g2 - 1] - noder[g1 - 1]) * xi
                self.setN0B0(N01, N02, B01, B02, xi, L, rr, cf, sf)
                tmp1 = rr * self.wt4[j] * dot(transpose(N01), N01)
                integ11 += tmp1
                tmp1 = rr * self.wt4[j] * dot(transpose(N01), N02)
                integ12 += tmp1
                tmp1 = rr * self.wt4[j] * dot(transpose(N02), N02)
                integ22 += tmp1
                tmp1 = rr * self.wt4[j] * dot(dot(transpose(B01), DD), B01)
                hnteg11 += tmp1
                tmp1 = rr * self.wt4[j] * dot(dot(transpose(B01), DD), B02)
                hnteg12 += tmp1
                tmp1 = rr * self.wt4[j] * dot(dot(transpose(B02), DD), B02)
                hnteg22 += tmp1
            #半分に割る(重み係数の分)
            integ11 /= 2
            integ12 /= 2
            integ22 /= 2
            hnteg11 /= 2
            hnteg12 /= 2
            hnteg22 /= 2
            ans11 = self.row * self.thick * 2 * math.pi * L * dot(dot(transpose(ramda), integ11), ramda)
            ans12 = self.row * self.thick * 2 * math.pi * L * dot(dot(transpose(ramda), integ12), ramda)
            ans22 = self.row * self.thick * 2 * math.pi * L * dot(dot(transpose(ramda), integ22), ramda)
            bns11 = 2 * math.pi * L * dot(dot(transpose(ramda), hnteg11), ramda)
            bns12 = 2 * math.pi * L * dot(dot(transpose(ramda), hnteg12), ramda)
            bns22 = 2 * math.pi * L * dot(dot(transpose(ramda), hnteg22), ramda)
            for j in xrange(3):
                for k in xrange(j, 3):
                    B[st1 + j][st1 + k] += bns11[j][k] - OMOM * ans11[j][k]
                    B[st2 + j][st2 + k] += bns22[j][k] - OMOM * ans22[j][k]
                    if k != j:
                        B[st1 + k][st1 + j] += bns11[j][k] - OMOM * ans11[j][k]
                        B[st2 + k][st2 + j] += bns22[j][k] - OMOM * ans22[j][k]
            for j in xrange(3):
                for k in xrange(3):
                    if st1 < st2:
                        B[st1 + j][st2 + k] += bns12[j][k] - OMOM * ans12[j][k]
                        if st1 + j != st2 + k:
                            B[st2 + k][st1 + j] += bns12[j][k] - OMOM * ans12[j][k]
                    else:
                        B[st2 + j][st1 + k] += bns12[k][j] - OMOM * ans12[k][j]
                        if st2 + j != st1 + k:
                            B[st1 + k][st2 + j] += bns12[k][j] - OMOM * ans12[k][j]
            if self.psz[i][0] == self.pointload:
                pass
            elif self.psz[i][0] == self.lineload:
                pass
            elif self.psz[i][0] == self.surfaceload:
                pass
            elif self.psz[i][0] == self.axisymln:
                integ12 = zeros((3, 2))
                hnteg12 = zeros((3, 2))
                if self.psz[i][1] == 1:
                    integ12[0][0] = integ12[1][1] = noder[g1 - 1] * 2 * math.pi
                elif self.psz[i][1] == 2:
                    hnteg12[0][0] = hnteg12[1][1] = noder[g2 - 1] * 2 * math.pi
                ans12 = dot(transpose(ramda), integ12)
                bns12 = dot(transpose(ramda), hnteg12)
                if speaker and i == 0:
                    for j in xrange(3):
                        B[st1+j][0] += ans12[j][0] * self.flux[i][0] * jomega * A + ans12[j][1] * self.flux[i][2] * jomega * A
                        B[st2+j][0] += bns12[j][0] * self.flux[i][0] * jomega * A + bns12[j][1] * self.flux[i][2] * jomega * A
                for j in xrange(3):
                    V[st1 + j] += ans12[j][0] * self.flux[i][0] + ans12[j][1] * self.flux[i][2]
                    V[st2 + j] += bns12[j][0] * self.flux[i][0] + bns12[j][1] * self.flux[i][2]
            elif self.psz[i][0] == self.axisymsf:
                pass
        return B, V
    
    def setN0B0(self, N01, N02, B01, B02, xi, L, rr, cf, sf):
        A1 = 1 - xi
        A2 = 1 - 2 * xi
        A3 = xi - xi**2
        A4 = 2 * xi - 3 * xi**2
        A5 = 1 - 4 * xi + 3 * xi**2
        A6 = xi**2 - xi**3
        A7 = 3 * xi**2 - 2 * xi**3
        A8 = 1 - 3 * xi**2 + 2 * xi**3
        A9 = xi - 2 * xi**2 + xi**3
        N01[0][0] = A1
        N01[1][1] = A8
        N01[1][2] = L * A9
        N02[0][0] = xi
        N02[1][1] = A7
        N02[1][2] = -L * A6
        B01[0][0] = -1 / L
        B01[1][0] = A1 * sf / rr
        B01[1][1] = A8 * cf / rr
        B01[1][2] = L * A9 * cf / rr
        B01[2][1] = 6 * A2 / L**2
        B01[2][2] = (4 - 6 * xi) / L
        B01[3][1] = 6 * A3 * sf / L / rr
        B01[3][2] = -A5 * sf / rr
        B02[0][0] = 1 / L
        B02[1][0] = xi * sf / rr
        B02[1][1] = A7 * cf / rr
        B02[1][2] = -L * A6 * cf / rr
        B02[2][1] = -6 * A2 / L**2
        B02[2][2] = (2 - 6 * xi) / L
        B02[3][1] = -6 * A3 * sf / L / rr
        B02[3][2] = A4 * sf / rr
    
    def setNnBn(self, Nn1, Nn2, Bn1, Bn2, xi, L, rr, cf, sf, moden):
        A1 = 1 - xi
        A2 = 1 - 2 * xi
        A3 = xi - xi**2
        A4 = 2 * xi - 3 * xi**2
        A5 = 1 - 4 * xi + 3 * xi**2
        A6 = xi**2 - xi**3
        A7 = 3 * xi**2 - 2 * xi**3
        A8 = 1 - 3 * xi**2 + 2 * xi**3
        A9 = xi - 2 * xi**2 + xi**3
        Nn1[0][0] = A1
        Nn1[1][1] = A1
        Nn1[2][2] = A8
        Nn1[2][3] = L * A9
        Nn2[0][0] = xi
        Nn2[1][1] = xi
        Nn2[2][2] = A7
        Nn2[2][3] = -L * A6
        Bn1[0][0] = -1 / L
        Bn1[1][0] = A1 * sf / rr
        Bn1[1][1] = moden * A1 / rr
        Bn1[1][2] = A8 * cf / rr
        Bn1[1][3] = L * A9 * cf / rr
        Bn1[2][0] = -moden * A1 / rr
        Bn1[2][1] = (-A1 * sf - rr / L) / rr
        Bn1[3][2] = 6 * A2 / L**2
        Bn1[3][3] = (4 - 6 * xi) / L
        Bn1[4][1] = moden * A1 * cf / rr**2
        Bn1[4][2] = (moden**2 * A8 / rr + 6 * A3 * sf / L) / rr
        Bn1[4][3] = (moden**2 * L * A9 / rr - A5 * sf) / rr
        Bn1[5][1] = -2 * (A1 * sf / rr + 1 / L) * cf / rr
        Bn1[5][2] = 2 * moden * (-6 * A3 / L - A8 * sf / rr) / rr
        Bn1[5][3] = 2 * moden * (A5 - L * A9 * sf / rr) / rr
        Bn2[0][0] = 1 / L
        Bn2[1][0] = xi * sf / rr
        Bn2[1][1] = moden * xi / rr
        Bn2[1][2] = A7 * cf / rr
        Bn2[1][3] = -L * A6 * cf / rr
        Bn2[2][0] = -moden * xi / rr
        Bn2[2][1] = (-xi * sf + rr / L) / rr
        Bn2[3][2] = -6 * A2 / L**2
        Bn2[3][3] = (2 - 6 * xi) / L
        Bn2[4][1] = moden * xi * cf / rr**2
        Bn2[4][2] = (moden**2 * A7 / rr - 6 * A3 * sf / L) / rr
        Bn2[4][3] = (-moden**2 * L * A6 / rr + A4 * sf) / rr
        Bn2[5][1] = -2 * (xi * sf / rr - 1 / L) * cf / rr
        Bn2[5][2] = 2 * moden * (6 * A3 / L - A7 * sf / rr) / rr
        Bn2[5][3] = 2 * moden * (-A4 + L * A6 * sf / rr) / rr
    
    def applybc0(self, B, V):
        ne = size(self.conec, 0)
        mnode3 = size(B, 0)
        mnode = mnode3 / 3
        for i in xrange(mnode): #loop for node
            cores = False
            for j in xrange(ne): #loop for element
                if self.psz[j][0] == self.ssupport: #check if element is simply supported
                    g1 = self.conec[j][0]
                    g2 = self.conec[j][1]
                    if (i == g1 - 1 and self.psz[j][1] == 1) or (i == g2 - 1 and self.psz[j][1] == 2):
                        st = i * 3
                        for k in xrange(st + 1):
                            B[k][st + 1] = B[st + 1][k] = 0
                        B[st + 1][st + 1] = 1
                        for k in xrange(st + 2, mnode3):
                            B[st + 1][k] = B[k][st + 1] = 0
                        V[st + 1] = 0
                        print ' node', i + 1, 'boundary condition is applied.'
                        cores = True
                elif self.psz[j][0] == self.psupport: #check if element is perfectly supported
                    g1 = self.conec[j][0]
                    g2 = self.conec[j][1]
                    if (i == g1 - 1 and self.psz[j][1] == 1) or (i == g2 - 1 and self.psz[j][1] == 2):
                        st = i * 3
                        for k in xrange(st):
                            B[k][st] = B[st][k] = B[k][st + 1] = B[st + 1][k] = 0
                        B[st][st] = 1
                        B[st][st + 1] = B[st + 1][st] = 0
                        B[st + 1][st + 1] = 1
                        for k in xrange(st + 2, mnode3):
                            B[st][k] = B[k][st] = B[st + 1][k] = B[k][st + 1] = 0
                        V[st] = V[st + 1] = 0
                        print ' node', i, 'boundary condition is applied.'
                        cores = True
    
    def calcuwb0(self, r, z, noder, nodez, V, eps):
        ne = size(self.conec, 0)
        un, wn, bn = 0, 0, 0
        for i in xrange(ne):
            g1 = self.conec[i][0]
            g2 = self.conec[i][1]
            st1 = (g1 - 1) * 3
            st2 = (g2 - 1) * 3
            r1 = noder[g1 - 1]
            z1 = nodez[g1 - 1]
            u1 = noder[g2 - 1] - r1
            v1 = nodez[g2 - 1] - z1
            L = sqrt(u1**2 + v1**2)
            if abs(u1) > eps and abs(v1) > eps:
                rxi = (r - r1) / u1
                zxi = (z - z1) / v1
                if rxi >= 0 and rxi <= 1 and abs(rxi - zxi) <= eps:
                    un, wn, bn = self.setuwb0(V, st1, st2, rxi, L)
                    break
            elif abs(u1) <= eps and abs(v1) > eps:
                zxi = (z - z1) / v1
                if zxi >= 0 and zxi <= 1 and abs(r - r1) <= eps:
                    un, wn, bn = self.setuwb0(V, st1, st2, zxi, L)
                    break
            elif abs(u1) > eps and abs(v1) <= eps:
                rxi = (r - r1) / u1
                if rxi >= 0 and rxi <= 1 and abs(z - z1) <= eps:
                    un, wn, bn = self.setuwb0(V, st1, st2, rxi, L)
                    break
        return un, wn, bn
    
    def setuwb0(self, V, st1, st2, xi, L):
        un = (1 - xi) * V[st1] + xi * V[st2]
        wn = (1 - 3 * xi**2 + 2 * xi**3) * V[st1 + 1] + L * (xi - 2 * xi**2 + xi**3) * V[st1 + 2]
        wn += (3 * xi**2 - 2 * xi**3) * V[st2 + 1] - L * (xi**2 - xi**3) * V[st2 + 2]
        bn = 6 / L * (xi - xi**2) * (V[st2 + 1] - V[st1 + 1]) + (1 - 4 * xi + 3 * xi**2) * V[st1 + 2] - (2 * xi - 3 * xi**2) * V[st2 + 2]
        return un, wn, bn
    
    def makeBV(self, BB, PP, noder, nodez, n, omega):
        pass
    
    def applybc(self, BB, PP):
        pass
    
    def calcuvwb(self, un, vn, wn, bn, r, z, noder, nodez, PP, eps):
        pass
    
    def spmkelement(self, divide, eps, noder, nodez):
        pass
    
    def spapplybc0(self, spBB, PP):
        pass
